%%
clc;
h1=[2,1,-3,5,4,5,-3,1,2];
h2=[2,1,-3,5,5,-3,1,2];
h3=[2,1,-3,5,0,-5,3,-1,-2];
h4=[2,1,-3,5,-5,3,-1,-2];
[H1,w1,type,tao]=amplres(h1);
type1=type;
tao1=tao;
subplot(221);
plot(w1/pi,H1);
grid;
title('线性相位Ⅰ型FIR滤波器幅度函数');
xlabel('\omega/\pi');
ylabel('H(\omega)');

[H2,w2,type,tao]=amplres(h2);
type2=type;
tao2=tao;
subplot(222);
plot(w2/pi,H2);
grid;
title('线性相位Ⅱ型FIR滤波器幅度函数');
xlabel('\omega/\pi');
ylabel('H(\omega)');

[H3,w3,type,tao]=amplres(h3);
type3=type;
tao3=tao;
subplot(223);
plot(w3/pi,H3);
grid;
title('非线性相位Ⅰ型FIR滤波器幅度函数');
xlabel('\omega/\pi');
ylabel('H(\omega)');

[H4,w4,type,tao]=amplres(h4);
type4=type;
tao4=tao;
subplot(224);
plot(w4/pi,H4);
grid;
title('非线性相位Ⅱ型FIR滤波器幅度函数');
xlabel('\omega/\pi');
ylabel('H(\omega)');

%%
wc=0.4*pi;
h=ideallp(wc,9);
%展示结果
[H,w]=freqz(h,1,512);
plot(w/pi,abs(H));
grid;
title('理想低通滤波器幅度特性');
xlabel('\omega/\pi');
ylabel('H(\omega)');

%%
%理想线性相位高通滤波器
wc=0.4*pi;
N=9;
hd=ideallp(pi, N)-ideallp(wc, N);
%展示结果
[H,w]=freqz(hd,1,512);
plot(w/pi,abs(H));
grid;
title('理想线性相位高通滤波器幅度特性');
xlabel('\omega/\pi');
ylabel('H(\omega)');

%%
%理想线性相位带通(带阻)滤波器，w1是通带上截止频率，w2是通带下截止频率
w1=0.3*pi;
w2=0.6*pi;
N=9;
%带通
hd1=ideallp(w2, N)-ideallp(w1, N);
%带阻
hd2=ideallp(pi, N)-hd1;
%展示结果
subplot(121);
[H1,w]=freqz(hd1,1,512);
plot(w/pi,abs(H1));
grid;
title('理想线性相位带通滤波器幅度特性');
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(122);
[H2,w]=freqz(hd2,1,512);
plot(w/pi,abs(H2));
grid;
title('理想线性相位带阻滤波器幅度特性');
xlabel('\omega/\pi');
ylabel('H(\omega)');

%%
%实现各类窗函数
N=33;
%矩形窗
w1=boxcar(N);
%三角窗
w2=triang(N);
%汉明窗
w3=hamming(N);
%汉宁窗
w4=hanning(N);
%布莱克曼窗
w5=blackman(N);
%凯泽窗
beta=0.5;
w6=kaiser(N,beta);
%展示结果
subplot(231);plot(w1);grid;title('矩形窗');
subplot(232);plot(w2);grid;title('三角窗');
subplot(233);plot(w3);grid;title('汉明窗');
subplot(234);plot(w4);grid;title('汉宁窗');
subplot(235);plot(w5);grid;title('布莱克曼窗');
subplot(236);plot(w6);grid;title('凯泽窗');

%%
clc;
wp=0.3*pi;
ws=0.45*pi;
deltaw=ws-wp;
N0=ceil(6.6*pi/deltaw);
N=N0+mod(N0+1,2);
n=0:N-1;
wd=(hamming(N))';
wc=(wp+ws)/2;
hd=ideallp(wc,N);
h=hd.*wd;
[db,mag,pha,grd,w]=freqz_m(h,1);
dw=2*pi/1000;
Rp=-min(db(1:wp/dw+1));
As=-max(db(ws/dw+1:501));
subplot(121);plot(n,wd);grid;title('汉明窗');xlabel('n');ylabel('w(n)');
axis([0,N,0,1]);

subplot(122);plot(w/pi,db);title('幅度响应（db）');xlabel('\omega/\pi');ylabel('20log|H(e^j^\omega)|(dB)');
axis([0,1,-80,5]);grid;

%%
% 设计一个线性相位FIR数字高通滤波器，技术指标为wp=0.7*pi,ws=0.5*pi,As=55dB，求h(n)并画出幅度响应（以dB表示）及相位响应曲线
wp = 0.7 * pi; 
ws = 0.5 * pi; 
As = 55; 

delta_w = wp - ws;
N = ceil((As - 7.95) / (2.285 * delta_w)) + 1;
if mod(N, 2) == 0
    N = N + 1;
end
wc = (wp + ws) / 2; 
h_lp = ideallp(wc, N);
h_hp = -h_lp;
h_hp((N-1)/2 + 1) = h_hp((N-1)/2 + 1) + 1;
window = hamming(N)';
h_hp = h_hp .* window;
[H, w] = freqz(h_hp, 1, 1024);

figure;
subplot(2, 1, 1);
plot(w/pi, 20*log10(abs(H)));
title('幅度响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('幅度 (dB)');
grid on;
subplot(2, 1, 2);
plot(w/pi, angle(H));
title('相位响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('相位 (弧度)');
grid on;

%%
%设计一个线性相位FIR数字高通滤波器，技术指标为wp1=0.4*pi,wp2=0.5*pi,ws1=0.2*pi,ws2=0.7*pi,As=75dB，求h(n)并画出幅度响应（以dB表示）及相位响应曲线
wp1 = 0.4 * pi; 
wp2 = 0.5 * pi; 
ws1 = 0.2 * pi; 
ws2 = 0.7 * pi; 
As = 75;    

delta_w = min(wp1 - ws1, ws2 - wp2);
N = ceil((As - 7.95) / (2.285 * delta_w)) + 1;
if mod(N, 2) == 0
    N = N + 1;
end
wc = (wp1 + wp2) / 2; 
h_lp = ideallp(wc, N);
h_hp = -h_lp;
h_hp((N-1)/2 + 1) = h_hp((N-1)/2 + 1) + 1;
window = hamming(N)';
h_hp = h_hp .* window;
[H, w] = freqz(h_hp, 1, 1024);

% 绘制幅度响应（以dB表示）
figure;
subplot(2, 1, 1);
plot(w/pi, 20*log10(abs(H)));
title('幅度响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('幅度 (dB)');
grid on;
% 绘制相位响应
subplot(2, 1, 2);
plot(w/pi, angle(H));
title('相位响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('相位 (弧度)');
grid on;

%%
% 设计一个线性相位FIR数字带阻滤波器。技术指标为wp1=0.35π，wp2=0.8π，wst1=0.5π，wst2=0.65π，As=80dB。求h(n)并画出幅度响应（以dB表示）及相位响应
wp1 = 0.35 * pi;  
wp2 = 0.8 * pi;  
ws1 = 0.5 * pi;  
ws2 = 0.65 * pi; 
As = 80;

delta_w = min(ws1 - wp1, wp2 - ws2);
N = ceil((As - 7.95) / (2.285 * delta_w)) + 1;
if mod(N, 2) == 0
    N = N + 1;
end
wc1 = (wp1 + ws1) / 2; 
wc2 = (wp2 + ws2) / 2; 
h_lp1 = ideallp(wc1, N); 
h_lp2 = ideallp(wc2, N); 
h_bp = h_lp1 + h_lp2;
h_bp((N-1)/2 + 1) = h_bp((N-1)/2 + 1) - 1;
window = hamming(N)';
h_bp = h_bp .* window;
[H, w] = freqz(h_bp, 1, 1024);

% 绘制幅度响应（以dB表示）
figure;
subplot(2, 1, 1);
plot(w/pi, 20*log10(abs(H)));
title('幅度响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('幅度 (dB)');
grid on;
% 绘制相位响应
subplot(2, 1, 2);
plot(w/pi, angle(H));
title('相位响应');
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('相位 (弧度)');
grid on;